• Modelling the spatial variability and the association with candidate driving factors in the reporting of violence against women (VAW)

• Variable of interest: counts of accesses to Anti-Violence Centers (AVCs) in Apulia per residence municipality of victims

• Observation period: 2021 – 2024; observation units: Apulian municipalities.

• Online version (recursive); more details here

• Offering shelter, legal assistance and psycological support to women victims of gender-based violence

• Maintaining systematic records of reported cases, as outlined in Italian Law no.119/2013

• Data collected by the Statistical Section and the Welfare Department of the Apulia Region, as part of a systematic monitoring program established under Regional Law no. 29/2014

• Distance (TEP) from the closest municipality hosting an AVC

• Components of the Municipality Frailty Index (MFI), 2021; source: ISTAT

• 256 municipalities, 4 years

• Standard Poisson linear regression: for \(i \in [1,256]\), and \(t \in \lbrace 2021, ..., 2024 \rbrace\) we have

\[ y_{it}\mid \eta_{it}, \theta \sim \mathrm{Poisson} \left( \mathrm{E}_{it} \, e^{\eta_{it}} \right) \]

• Where

  • \(y_{it}\): AVC accesses in municipality \(i\) and year \(t\);

  • \(\mathrm{E}_{it}\): female population aged \(>14\) years (offset)

  • \(\eta_{it}\): linear predictor; \(\theta\): model hyperparameters

• Inclusion of a multivariate latent spatial effect in \(\eta\) to account for territorial variability and overdispersion

\[\eta_{it} := X_{it}^\top \beta + Z_{it}\] - \(X\): Covariates; \(\beta\): covariate effects

• \(Z_{it}\): multivariate latent spatial effect of size \(n \times k\), where \(n=256\) and \(k=4 \quad \longrightarrow\) modelled via the flexible M-Models (Botella-Rocamora, Martinez-Beneito, and Banerjee 2015)

• Competing models: ICAR, PCAR, LCAR and BYM

• Simplest formulation among spatial autoregressive models

\[ Z \mid \Sigma \sim \mathcal{N}_{kn} \left( 0, \Sigma \otimes (D - W)^{+} \right) \] - \(\Sigma\): \(4 \times 4\) between-fields variance-covariance matrix

• \(D\): Graph degree matrix;

• \(W\): Graph neighbourhood matrix

\(+\) Easy to fit and interpret

\(-\) implicitly assumes spatial autocorrelation equal to 1 \(\rightarrow\) sometimes too smooth to adequately fit data

• Based on the decomposition of the variance-covariance matrix as \(\Sigma := M^\top M\)

• \(M\): full rank matrix, of size \(k \times k\) herein

• consider \(Z\) as a combination of \(k\) independent fields \(Z_j^{(0)}\) for \(1 \leq j \leq k\), i.e. \(Z := \left(Z_1^{(0)}, \dots, Z_k^{(0)}\right) M\)

• \(Z^{(0)} \sim \mathcal{N}_n \left( 0, \Omega_j \right)\) for a \(n \times n\) s.p.d. matrix \(\Omega_j\) depending:

  • On the graph structure
  • On an array of model-specific parameters

Joint prior of \(Z\):

\[ Z \sim \mathcal{N}_{nk} \left(0, (M^\top\otimes I_n)\, \Omega \, (M \otimes I_n) \right) \] Where \(\Omega\) is block-diagonal: \[ \Omega = \begin{pmatrix} \Omega_1 & 0 & \dots & 0 \\ 0 & \Omega_2 & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & \Omega_k \end{pmatrix} \]

• Three standard models for \(\Omega\) compared in next slides

• Generalises the ICAR allowing for a proper prior density with nonsingular precision (Gelfand and Vounatsou 2003)

• Introduces a spatial correlation parameter \(A := \mathrm{diag}(\alpha_1, ... \alpha_k)\), \(\alpha_j \in [0,1] \, \forall j\)

\[ \Omega = (I_k \otimes D - A \otimes W)^{-1} \]

• Limit cases:
  • \(\alpha_t = 0\): independent model (\(\neq\) IID)
  • \(\alpha_t = 1\): ICAR

• Also generalises the ICAR allowing for a proper density

• Introduces a precision mixing parameter \(\Lambda := \mathrm{diag}(\lambda_1, ... \lambda_k)\), yet with \(\lambda_j \in [0,1] \, \forall j\) (Leroux, Lei, and Breslow 2000)

\[ \Omega_j = \left(\Lambda \otimes L + (I_k - \Lambda) \right) \] - \(L\) = graph Laplacian matrix, i.e. \(L = D-W\)

• Limit cases:
  • \(\lambda_t = 0\): IID model
    
  • \(\lambda_t = 1\): ICAR

• Actually a reparametrisation of the original BYM (Besag, York, and Mollié 1991), proposed by (Riebler et al. 2016)

\[ Z = U \Phi^\frac{1}{2} M + V (I_k-\Phi^\frac{1}{2}) M \] - \(\Phi := \mathrm{diag}(\phi_1, \dots, \phi_k)\): variance mixing parameter, \(\phi_j \in [0, 1]\, \forall j\)

• \(U \sim \mathcal{N}_{kn} \left(0, I_k \otimes L_S^{+} \right)\): independent ICAR (limit case for \(\Phi = I_k\)) \(\rightarrow\) spatial dependence, where \(L_S\): scaled Laplacian (Riebler et al. 2016)

• \(V \sim \mathcal{N}_{kn} (0, I_{kn})\): IID field (limit case for \(\Phi = 0\)) \(\rightarrow\) noise

• The parameters \(A\), \(\Lambda\) and \(\Phi\) are usually assigned a Uniform prior

• A more principle-based approach is using PC-priors

  • This was already done by Rieber et Al. (2016) for the univariate BYM

• Question: can we apply this approach in the multivariate setting, and also for PCAR and LCAR models?

• Caveat: deriving analytically a PC-prior for multivariate parameters is hardly possible (Simpson et al. 2017), (Battagliese 2020) \(\longrightarrow\) workaround: Sequential derivation

Penalise the departure of a model depending on a parameter \(\theta\) from a simpler (base) one (Simpson et al. 2017):

• Simple model: here, an independent, unstructured one, with \(\theta = 0\)

• Departure: a measure of how the more flexible model, say \(\mathcal{M}_F\) diverges from the simple one, say \(\mathcal{M}_B \rightarrow\) distance function: \(d (\mathcal{M}_F||\mathcal{M}_B) := \sqrt{2 \mathrm{KLD}(\mathcal{M}_F||\mathcal{M}_B)}\)

• Penalise: shrink the p.d.f. of \(d(\mathcal{M}_F||\mathcal{M}_B)\) at a constant rate \(\rightarrow\) exponential prior on \(d\)

Fix two parameters \(\theta_0, \mathcal{a}:\mathrm{Prob} \lbrace \theta \leq \theta_0 \rbrace = \mathcal{a}\)

• Let \(\theta_j\) denote either \(\alpha_j\), \(\lambda_j\), or \(\phi_j\) for \(j = 1, \dots, k\)

• At each step \(j\), we test the flexible model with \(\theta_j > 0\), against the base model with \(\theta_j = 0\), all other parameters being unchanged

• Then, the joint prior simply comes from the chain rule:

\[\pi(\theta) = \pi (\theta_1) \pi(\theta_2 \mid \theta_1) \dots \pi(\theta_k \mid \theta_1, \dots, \theta_{k-1})\]

• Provided that

\[ \begin{align*} \pi(\theta_j \mid \theta_1, \dots, \theta_{j-1}) = \\ =\pi(\theta_j \mid \theta_1, \dots, \theta_{j-1}, \theta_{j+1}=\dots=\theta_k=0) \end{align*} \]

• PC-priors available in closed formula

• \(\pi(\theta)\) independent on \(\Sigma\)

• \(\pi(\theta_j \mid \theta_1, \dots, \theta_{j-1}) = \pi(\theta_j)\), i.e., at each step \(j\), the prior of \(\theta_j\) is independent on other parameters

\[ \therefore \pi(\theta) = \prod_{j=1}^k \pi(\theta_j) \]

• Only needed Laplacian eigenvalues (LCAR, BYM) or eigenvalues of \(D^{-1}W\) (PCAR), both deterministic \(\rightarrow\) computed once

Models are compared by the Watanabe-Akaike Information Criterion (WAIC) (Gelman, Hwang, and Vehtari 2014). Free parameters are reported as well.

• Covariate effects: independent weekly informative Gaussians \(\mathcal{N}(0, 10^3)\)

• \(\Sigma\): Wishart prior, with identity scale parameter and 8 degrees of freedom

• \(\Phi\): PC-prior, with \(\mathrm{Prob} \lbrace \phi_j \leq 0.6 \rbrace = 0.9\) for all \(j\)

• Model estimation is carried out entirely within R-INLA (Rue, Martino, and Chopin 2009) using the Variational Bayes (Van Niekerk et al. 2023) correction to the posterior mean of latent variables

\(y =\) VAW occurrence \(\times\) incidents reporting rate

• Twofold interpretation of covariate effects on AVC accesses

• TEP: either relative ease in reporting violence in proximity of AVCs, or appropriateness of AVC placement

• ELL: lower reliance on AVCs in poor education contexts

• ER: lower violence occurrence in more developed areas

• Posterior distribution of the mixing parameter (\(\Phi\))

• Percentage of variation explained by spatial structure

• Posterior of the variance-covariance matrix (\(\Sigma\))

Posterior values of \(\mathbb{E}[Z \mid y]\)